import geopandas as gpd
import pandas as pd  # 新增pandas导入
import numpy as np
import os
import gc
import math
from shapely.geometry import Polygon, MultiPolygon
from shapely.validation import make_valid
import re
from datetime import datetime

def export_forest_model_clean(geojson_path, output_obj_path, scale=100000, thickness=2.0):
    """
    适配「鹤北耕地」GeoJSON的模型导出：
    1. 匹配GeoJSON中的「林班号/小班号/林场」字段
    2. 优化MultiPolygon处理，避免几何破碎
    3. 移除无效面/碎面，保留完整3D厚度
    """
    print(f"🌲 开始导出适配'鹤北耕地'的优化模型...")
    print(f"📁 输入GeoJSON: {geojson_path}")
    print(f"📁 输出OBJ文件: {output_obj_path}")
    print(f"📏 厚度设置: {thickness} 单位")
    print(f"✨ 优化目标：适配字段+移除无效面+紧密贴合")
    print(f"⏰ 开始时间: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")

    # 检查输入文件
    if not os.path.exists(geojson_path):
        print(f"❌ 错误：输入文件不存在 - {geojson_path}")
        return False

    try:
        # 1. 读取并预处理数据（适配鹤北耕地GeoJSON字段）
        gdf = gpd.read_file(geojson_path)
        print(f"✅ 成功读取GeoJSON：共包含 {len(gdf)} 个地块")

        # 查看字段名，确保匹配（适配「林班号/小班号/林场」字段）
        required_fields = ["林班号", "小班号", "林场"]
        for field in required_fields:
            if field not in gdf.columns:
                print(f"⚠️  警告：GeoJSON中未找到'{field}'字段，将使用默认值")

        # 数据清洗与标准化
        valid_plots = []
        for idx, row in gdf.iterrows():
            # 适配GeoJSON字段：林班号（数字）、小班号（数字）、林场（字符串）
            林场 = str(row.get("林场", "未知林场")).strip()
            # 修复pd未定义错误，处理空值
            林班_val = row.get("林班号", 0)
            林班 = str(int(林班_val) if not pd.isna(林班_val) else "0000").strip().zfill(4)

            小班_val = row.get("小班号", 0)
            小班 = str(int(小班_val) if not pd.isna(小班_val) else "00000").strip().zfill(5)

            几何 = row.geometry

            # 清理林场名称特殊字符（适配Windows和Unity）
            林场 = re.sub(r'[\\/*?:"<>|]', '', 林场).replace(' ', '_')
            完整名称 = f"{林场}_{林班}_{小班}"

            # 修复无效几何（重点处理MultiPolygon）
            if not 几何.is_valid:
                几何 = make_valid(几何)
                print(f"⚠️  修复无效几何：{完整名称}（原类型：{type(几何).__name__}）")

            # 过滤空几何和微小地块（面积阈值适配耕地尺度）
            if not 几何.is_empty:
                # 计算地块面积（基于CRS84经纬度的近似面积，单位：平方米）
                geom_area = 几何.area * 111000 * 111000 * math.cos(math.radians(几何.centroid.y))
                if geom_area > 10:  # 过滤小于10平方米的微小地块
                    valid_plots.append({
                        "林场": 林场,
                        "林班": 林班,
                        "小班": 小班,
                        "完整名称": 完整名称,
                        "几何": 几何,
                        "面积": geom_area
                    })
                else:
                    print(f"⚠️  跳过微小地块（{geom_area:.1f}㎡）：{完整名称}")
            else:
                print(f"⚠️  跳过空几何：{完整名称}")

        if not valid_plots:
            print(f"❌ 无有效地块数据")
            return False
        print(f"✅ 有效地块数量：{len(valid_plots)}（面积≥10㎡）")

        # 2. 建立层级结构（林场→林班→小班）
        forest_hierarchy = {}
        for plot in valid_plots:
            林场 = plot["林场"]
            林班 = plot["林班"]

            if 林场 not in forest_hierarchy:
                forest_hierarchy[林场] = {}
            if 林班 not in forest_hierarchy[林场]:
                forest_hierarchy[林场][林班] = []

            forest_hierarchy[林场][林班].append(plot)

        print(f"✅ 层级结构：{len(forest_hierarchy)} 个林场 → {sum(len(lb) for lb in forest_hierarchy.values())} 个林班")

        # 3. 计算全局偏移（适配CRS84经纬度，避免大坐标误差）
        all_coords = []
        for plot in valid_plots:
            geom = plot["几何"]
            processed_geoms = process_geometry(geom)  # 专门处理MultiPolygon
            for g in processed_geoms:
                # 提取经纬度坐标（忽略z值，GeoJSON中z=0.0可忽略）
                exterior_coords = [(coord[0], coord[1]) for coord in g.exterior.coords]
                all_coords.extend(exterior_coords)
                for interior in g.interiors:
                    interior_coords = [(coord[0], coord[1]) for coord in interior.coords]
                    all_coords.extend(interior_coords)

        coords_array = np.array(all_coords, dtype=np.float64)
        min_lon, min_lat = np.min(coords_array, axis=0)
        max_lon, max_lat = np.max(coords_array, axis=0)
        global_offset = np.array([(min_lon + max_lon) / 2, (min_lat + max_lat) / 2], dtype=np.float64)
        print(f"📊 全局地理范围（CRS84经纬度）：\n   经度：{min_lon:.8f} ~ {max_lon:.8f}\n   纬度：{min_lat:.8f} ~ {max_lat:.8f}")

        # 4. 创建坐标映射表（确保相邻地块共享顶点，消除间隙）
        print(f"🔍 创建坐标映射表...")
        coord_map = create_coord_map(all_coords, precision=8)  # 经纬度8位精度足够
        print(f"✅ 坐标映射表创建完成：{len(coord_map)} 个唯一坐标点")

        # 5. 模型参数配置（适配耕地模型）
        COORDINATE_PRECISION = 8    # 匹配经纬度精度
        MIN_FACE_AREA = 0.1         # 过滤0.1㎡以下的碎面

        # 6. 生成OBJ数据
        obj_data = []
        total_vertices = 0
        total_faces = 0
        removed_faces = 0

        # 6.1 文件头（标注适配信息）
        obj_data.append(f"# 鹤北耕地模型（优化版）")
        obj_data.append(f"# 适配CRS84经纬度 | 移除无效面 | 3D厚度 {thickness}")
        obj_data.append(f"# 生成时间：{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")
        obj_data.append(f"# 统计：{len(forest_hierarchy)} 林场 | {sum(len(lb) for lb in forest_hierarchy.values())} 林班 | {len(valid_plots)} 地块")
        obj_data.append("")

        # 6.2 写入顶点（顶部z=0 + 底部z=-thickness，垂坠厚度）
        print(f"📝 写入顶点数据...")
        sorted_coords = sorted(coord_map.items(), key=lambda x: x[1])
        vertices = []  # 存储顶点坐标，用于面积计算
        for (lon, lat), coord_idx in sorted_coords:
            # 经纬度→平面坐标转换（偏移+缩放，消除大坐标）
            x = (lon - global_offset[0]) * scale
            y = (lat - global_offset[1]) * scale

            # 顶部顶点（地表，z=0）
            top_v = (x, y, 0.0)
            obj_data.append(f"v {x:.{COORDINATE_PRECISION}f} {y:.{COORDINATE_PRECISION}f} 0.0")
            vertices.append(top_v)

            # 底部顶点（垂坠，z=-thickness）
            bottom_v = (x, y, -thickness)
            obj_data.append(f"v {x:.{COORDINATE_PRECISION}f} {y:.{COORDINATE_PRECISION}f} {-thickness}")
            vertices.append(bottom_v)

        total_vertices = len(vertices)
        obj_data.append(f"\n# 顶点总数：{total_vertices}（每个坐标点2个顶点）")
        obj_data.append("")

        # 6.3 生成优化后的面（移除无效面，处理MultiPolygon）
        print(f"🔗 生成优化面（移除碎面/重复面）...")
        for 林场名, 林班字典 in forest_hierarchy.items():
            # 林场层级（OBJ对象标签 o）
            obj_data.append(f"o {林场名}")
            obj_data.append(f"usemtl Material_{林场名}")
            obj_data.append(f"# 林场：{林场名}（{len(林班字典)}个林班）")
            obj_data.append("")

            for 林班名, 地块列表 in 林班字典.items():
                # 林班层级（OBJ组标签 g）
                林班组名 = f"{林场名}_{林班名}"
                obj_data.append(f"g {林班组名}")
                obj_data.append(f"usemtl Material_{林班组名}")
                obj_data.append(f"# 林班：{林班组名}（{len(地块列表)}个地块）")
                obj_data.append("")

                for 地块 in 地块列表:
                    # 地块层级（OBJ子组标签 g）
                    地块名 = 地块["完整名称"]
                    obj_data.append(f"g {地块名}")
                    obj_data.append(f"usemtl Material_{地块名}")
                    obj_data.append(f"# 地块：{地块名}（面积：{地块['面积']:.1f}㎡）")

                    # 处理MultiPolygon拆分后的单个Polygon
                    processed_geoms = process_geometry(地块["几何"])
                    for geom_idx, g in enumerate(processed_geoms):
                        # 生成面并过滤无效面
                        faces, removed = generate_optimized_faces(
                            g, coord_map, global_offset, scale, thickness,
                            vertices=vertices, min_area=MIN_FACE_AREA, precision=COORDINATE_PRECISION
                        )
                        removed_faces += removed

                        if not faces:
                            print(f"⚠️  地块{地块名}_part{geom_idx}无有效面，跳过")
                            continue

                        # 写入面到OBJ
                        for face in faces:
                            obj_data.append(f"f {face[0]} {face[1]} {face[2]}")
                        total_faces += len(faces)

                    obj_data.append("")  # 地块间分隔
                obj_data.append("")  # 林班间分隔
            obj_data.append("")  # 林场间分隔

        # 7. 写入OBJ文件
        os.makedirs(os.path.dirname(output_obj_path), exist_ok=True)
        with open(output_obj_path, 'w', encoding='utf-8') as f:
            f.write('\n'.join(obj_data))

        # 8. 输出最终统计
        print(f"\n🎉 鹤北耕地模型导出完成！")
        print(f"📊 模型统计：")
        print(f"   - 层级结构：{len(forest_hierarchy)} 林场 → {sum(len(lb) for lb in forest_hierarchy.values())} 林班 → {len(valid_plots)} 地块")
        print(f"   - 几何数据：{total_vertices:,} 个顶点 | {total_faces:,} 个有效面 | 移除 {removed_faces} 个无效面")
        print(f"   - 模型参数：厚度 {thickness} 单位 | 坐标精度 {COORDINATE_PRECISION} 位 | 最小面面积 {MIN_FACE_AREA} ㎡")
        print(f"   - 文件路径：{output_obj_path}")
        print(f"⏰ 完成时间: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")

        gc.collect()
        return True

    except Exception as e:
        print(f"❌ 导出错误: {str(e)}")
        import traceback
        traceback.print_exc()
        return False

def create_coord_map(coords, precision=8):
    """创建经纬度坐标映射表，确保相邻地块共享顶点"""
    coord_map = {}
    idx = 0
    for (lon, lat) in coords:
        # 经纬度四舍五入到指定精度，消除浮点误差
        lon_rounded = round(lon, precision)
        lat_rounded = round(lat, precision)
        key = (lon_rounded, lat_rounded)

        if key not in coord_map:
            coord_map[key] = idx
            idx += 1
    return coord_map

def process_geometry(geom):
    """专门处理MultiPolygon，拆分为单个Polygon，确保几何完整性"""
    processed = []
    if isinstance(geom, Polygon):
        if not geom.is_empty:
            # 确保Polygon闭合（处理GeoJSON中可能的未闭合问题）
            if not np.allclose(geom.exterior.coords[0][:2], geom.exterior.coords[-1][:2], atol=1e-8):
                exterior_coords = list(geom.exterior.coords) + [geom.exterior.coords[0]]
                # 保留内部孔洞
                interiors = [list(interior.coords) for interior in geom.interiors]
                closed_poly = Polygon(exterior_coords, interiors)
                processed.append(closed_poly)
            else:
                processed.append(geom)
    elif isinstance(geom, MultiPolygon):
        # 拆分MultiPolygon为单个Polygon
        for poly in geom.geoms:
            if not poly.is_empty:
                if not np.allclose(poly.exterior.coords[0][:2], poly.exterior.coords[-1][:2], atol=1e-8):
                    exterior_coords = list(poly.exterior.coords) + [poly.exterior.coords[0]]
                    interiors = [list(interior.coords) for interior in poly.interiors]
                    closed_poly = Polygon(exterior_coords, interiors)
                    processed.append(closed_poly)
                else:
                    processed.append(poly)
    return processed

def generate_optimized_faces(geom, coord_map, global_offset, scale, thickness, vertices, min_area=0.1, precision=8):
    """生成优化面：过滤碎面/重复面，处理经纬度坐标"""
    faces = []
    removed = 0
    face_set = set()  # 去重：存储排序后的面索引

    # 提取当前Polygon的经纬度坐标（忽略z值）
    exterior_coords = [(coord[0], coord[1]) for coord in geom.exterior.coords]
    # 映射为全局坐标索引（顶部+底部）
    top_indices = []  # 顶部顶点索引（z=0）
    bottom_indices = []  # 底部顶点索引（z=-thickness）
    for (lon, lat) in exterior_coords:
        lon_rounded = round(lon, precision)
        lat_rounded = round(lat, precision)
        key = (lon_rounded, lat_rounded)

        if key not in coord_map:
            continue  # 理论上不会出现，坐标已提前收集

        coord_idx = coord_map[key]
        # 顶部顶点索引 = coord_idx * 2 + 1（OBJ索引从1开始）
        top_idx = coord_idx * 2 + 1
        # 底部顶点索引 = coord_idx * 2 + 2
        bottom_idx = coord_idx * 2 + 2

        top_indices.append(top_idx)
        bottom_indices.append(bottom_idx)

    # 1. 生成顶部面（地表，过滤碎面）
    if len(top_indices) >= 3:
        # 三角化：从第一个顶点发散
        for i in range(1, len(top_indices) - 1):
            face = (top_indices[0], top_indices[i], top_indices[i+1])
            # 去重：排序后判断
            face_sorted = tuple(sorted(face))
            if face_sorted in face_set:
                removed += 1
                continue
            # 计算面面积，过滤过小的面
            face_area = calculate_face_area(face, vertices)
            if face_area < min_area:
                removed += 1
                continue
            # 有效面：添加到列表
            faces.append(face)
            face_set.add(face_sorted)

        # 2. 生成底部面（垂坠面，反向三角化避免背面剔除）
        for i in range(len(bottom_indices) - 2, 0, -1):
            face = (bottom_indices[0], bottom_indices[i], bottom_indices[i+1])
            face_sorted = tuple(sorted(face))
            if face_sorted in face_set:
                removed += 1
                continue
            face_area = calculate_face_area(face, vertices)
            if face_area < min_area:
                removed += 1
                continue
            faces.append(face)
            face_set.add(face_sorted)

        # 3. 生成侧面（连接顶部和底部，形成厚度）
        for i in range(len(top_indices) - 1):
            v1_top = top_indices[i]
            v2_top = top_indices[i+1]
            v2_bottom = bottom_indices[i+1]
            v1_bottom = bottom_indices[i]

            # 侧面第一个三角面
            face1 = (v1_top, v2_top, v2_bottom)
            face1_sorted = tuple(sorted(face1))
            if face1_sorted not in face_set and calculate_face_area(face1, vertices) >= min_area:
                faces.append(face1)
                face_set.add(face1_sorted)
            else:
                removed += 1

            # 侧面第二个三角面
            face2 = (v1_top, v2_bottom, v1_bottom)
            face2_sorted = tuple(sorted(face2))
            if face2_sorted not in face_set and calculate_face_area(face2, vertices) >= min_area:
                faces.append(face2)
                face_set.add(face2_sorted)
            else:
                removed += 1

    # 处理内部孔洞（与外部轮廓反向）
    for interior in geom.interiors:
        interior_coords = [(coord[0], coord[1]) for coord in interior.coords]
        hole_top = []
        hole_bottom = []
        for (lon, lat) in interior_coords:
            lon_rounded = round(lon, precision)
            lat_rounded = round(lat, precision)
            key = (lon_rounded, lat_rounded)
            if key not in coord_map:
                continue
            coord_idx = coord_map[key]
            hole_top.append(coord_idx * 2 + 1)
            hole_bottom.append(coord_idx * 2 + 2)

        if len(hole_top) >= 3:
            # 孔洞顶部面（反向三角化）
            for i in range(len(hole_top) - 2, 0, -1):
                face = (hole_top[0], hole_top[i], hole_top[i+1])
                face_sorted = tuple(sorted(face))
                if face_sorted in face_set or calculate_face_area(face, vertices) < min_area:
                    removed += 1
                    continue
                faces.append(face)
                face_set.add(face_sorted)

            # 孔洞底部面（正向三角化）
            for i in range(1, len(hole_bottom) - 1):
                face = (hole_bottom[0], hole_bottom[i], hole_bottom[i+1])
                face_sorted = tuple(sorted(face))
                if face_sorted in face_set or calculate_face_area(face, vertices) < min_area:
                    removed += 1
                    continue
                faces.append(face)
                face_set.add(face_sorted)

            # 孔洞侧面
            for i in range(len(hole_top) - 1):
                v1_top = hole_top[i]
                v2_top = hole_top[i+1]
                v2_bottom = hole_bottom[i+1]
                v1_bottom = hole_bottom[i]

                face1 = (v2_top, v1_top, v2_bottom)
                face1_sorted = tuple(sorted(face1))
                if face1_sorted not in face_set and calculate_face_area(face1, vertices) >= min_area:
                    faces.append(face1)
                    face_set.add(face1_sorted)
                else:
                    removed += 1

                face2 = (v2_bottom, v1_top, v1_bottom)
                face2_sorted = tuple(sorted(face2))
                if face2_sorted not in face_set and calculate_face_area(face2, vertices) >= min_area:
                    faces.append(face2)
                    face_set.add(face2_sorted)
                else:
                    removed += 1

    return faces, removed

def calculate_face_area(face, vertices):
    """计算三角面的实际面积（平方米）"""
    try:
        # 从顶点索引获取坐标（OBJ索引-1 = 列表索引）
        v1 = vertices[face[0] - 1]
        v2 = vertices[face[1] - 1]
        v3 = vertices[face[2] - 1]

        # 提取x,y坐标（平面坐标）
        x1, y1, _ = v1
        x2, y2, _ = v2
        x3, y3, _ = v3

        # 计算三角形面积（ shoelace公式 ）
        area = 0.5 * abs((x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1))
        return area
    except Exception as e:
        # 异常时返回0，确保被过滤
        return 0

# 执行入口（适配鹤北耕地GeoJSON路径）
if __name__ == "__main__":
    # 配置路径（替换为你的实际路径）
    input_geojson = r"C:\Users\ma\Desktop\鹤北geojson\鹤北耕地.geojson"
    output_obj = r"C:\Users\ma\Desktop\鹤北geojson\鹤北耕地模型_优化版.obj"

    # 执行导出（可调整厚度）
    success = export_forest_model_clean(
        geojson_path=input_geojson,
        output_obj_path=output_obj,
        scale=100000,  # 经纬度缩放系数，无需修改
        thickness=3.0  # 厚度调整为3.0，更明显
    )

    if success:
        print("\n🎉 鹤北耕地模型导出完成！模型无破碎感，相邻地块紧密贴合")
    else:
        print("\n❌ 导出失败，请查看错误信息")
